subroutine ReadMeshOpenFOAM()
!use Acoustics, only: AcousticSurf
use Data_cfoam
use Data_Mesh
!use Data_Transducer !, only: defCellsTransducer
implicit none
    include 'mpif.h'

    integer :: iconvname
    
    integer :: i,j,k, m,l, jj,j1
    integer :: isdispls, jstart, jpropeller
    integer :: netr0
    integer :: ier
    integer :: ndeltap
    
    real(8) :: xg(3), r
    real(8) :: vol ! volume
    real(8) :: Vmax, Vmin ! max and min volume
    real(8) :: Dbmin,Dfmin ! corrector to current and next face
    real(8) :: ds1

    integer :: iFace, PointNum, iPoint, FaceNum
    integer :: iProc,statusp(mpi_status_size)       
	! Read Boundary Mesh / Boundary Mesh has FOAM's format
	call inputMain(processor(1:ie)//'constant/polyMesh/owner',nPoints,nCells,nFaces,nInternalFaces)


    call MPI_ALLREDUCE(nPoints,nPointsGlobal,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
    call MPI_ALLREDUCE(nCells,nCellsGlobal,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
    call MPI_ALLREDUCE(nFaces,nFacesGlobal,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
    call MPI_ALLREDUCE(nInternalFaces,nInternalFacesGlobal,1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,ierr)
    
    if(rank.eq.0)then
	    write(*,*) 'nPoints (with proc duplicates):', nPointsGlobal
	    write(*,*) 'nCells:', nCellsGlobal
	    write(*,*) 'nFaces (with proc duplicates):', nFacesGlobal
	    write(*,*) 'nInternalFaces:', nInternalFacesGlobal
    end if

	
	call inputBoundary(processor(1:ie)//'constant/polyMesh/boundary',nBoundaryFields)
	allocate (boundaryMeshSize(nBoundaryFields) &
                ,boundaryMeshStart(nBoundaryFields) &
                ,boundaryMeshName(nBoundaryFields) &
                ,boundaryMeshType(nBoundaryFields) &
                ,BoundaryFieldFactor(nBoundaryFields) &
                ,nodeBoundaryMeshSize(nBoundaryFields) &
                ,nodeBoundaryMeshStart(nBoundaryFields) &
                ) 
    ! Read boundary Mesh
	call inputBoundary1(processor(1:ie)//'constant/polyMesh/boundary',nBoundaryFields &
                ,boundaryMeshSize                                        &
                ,boundaryMeshStart                                       &
                ,boundaryMeshName,boundaryMeshType)

    BoundaryFieldFactor=2.0d0 ! not used
    
    ! Set boundary between processors' fields
    do j=1,nBoundaryFields
        if (boundaryMeshName(j)(1:4).eq.'proc') then
	        BoundaryFieldFactor(j)=1.0d0
!	        write(*,*) 'MeshType=', boundaryMeshType(j), rank
!	        write(*,*) 'MeshName=', boundaryMeshName(j), rank
!	        write(*,*) 'MeshSize=', boundaryMeshSize(j),boundaryMeshStart(j), rank
!	        do j1=1,18
!	            if (boundaryMeshName(j)(j1:j1+1).eq.'to') &
!	                read(boundaryMeshName(j)(j1+2:20),*) BoundaryProcessorNumber(j)  ! Semiletov: It's a rude mistake, BoundaryProcessorNumber(nump+1) but nump+1 may be less nBoundaryFields
!!	            end if
!	        end do
	    end if
    !    if (boundaryMeshName(j)(1:4).eq.'INLT') BoundaryFieldFactor(j)=1.0d0
        if (boundaryMeshName(j)(1:4).eq.'OUTL') BoundaryFieldFactor(j)=1.0d0
    end do
    ! End of Setting boundary between processors' fields

    ! Read Interior Mesh / Interior Mesh has FOAM's format
	allocate (meshFaceOwner(nFaces))
	call inputMain1(processor(1:ie)//'constant/polyMesh/owner',nFaces,meshFaceOwner)
	
	allocate (meshFaceNeighbour(nInternalFaces))
	call inputMain1(processor(1:ie)//'constant/polyMesh/neighbour',nInternalFaces,meshFaceNeighbour)
	allocate (meshFaces(4,nFaces))	
	call inputMain2(processor(1:ie)//'constant/polyMesh/faces',nFaces,meshFaces)
	allocate (meshPoints(3,nPoints),pointProcAddressing(nPoints))
	! End. Read Interior Mesh / Interior Mesh has FOAM's format
	! MPI
!	allocate (pointLocalAddressing(nPoints))
	if (nump.gt.1) then
!	    open(11,file='processor'//char48(is:4)//'/constant/polyMesh/pointProcAddressing')
!        write(*,*) processor(1:ie)//'constant/polyMesh/pointProcAddressing'
        open(11,file=processor(1:ie)//'constant/polyMesh/pointProcAddressing',iostat=ierr)
        if(ierr.ne.0) then
            write(*,*) '@@ Fatal error: ', processor(1:ie)//'constant/polyMesh/pointProcAddressing', &
                       ' open error!'
            stop '@@ Fatal error'
        end if

!	    open(11,file=processor//'constant/polyMesh/pointProcAddressing')
	            ! 'processor'//char48(is:4)//'/constant/polyMesh/pointProcAddressing')	    
	    do i=1,18
	        read(11,*)
	    end do
	    read (11,*) 
	    read(11,*)
	    do i=1,nPoints
	        read(11,*) pointProcAddressing(i)
	    end do
	    close(11)
	else
	    do i=1,nPoints
	        pointProcAddressing(i)=i
	    end do
	end if

    call CalcNodeInterface()

	! MPI + rotating coor system (propeller)
	call inputMain3(processor(1:ie)//'constant/polyMesh/points',nPoints,meshPoints)
	
	do i = 1,nPoints
	    do j = 1,3
    	    meshPoints(j,i) = meshFactor(j) * meshPoints(j,i)
    	end do
	end do
	
!	do i = 1,nPoints
!    	meshPoints(1,i) = meshPoints(1,i) + 0.5*meshPoints(2,i) - 25d0
!	end do
	
	
	
!    do i = 1,nodeInterface
!        write(2,'(a5,i5,a3,i8,a6,i8,a7,3f10.5)') 'rank=', rank, ' i=', i, ' Node=', NodeBoundaryMeshLabel(i), ' Coord=', meshPoints(1:3,NodeBoundaryMeshLabel(i))
!    end do
	
	isdispls=-1
	jstart=-1
	jpropeller=0
	do k=1,nBoundaryFields
!        if (boundaryMeshName(k)(1:4).eq.'WINT') jpropeller=k
	    if (boundaryMeshName(k)(1:4).eq.'proc') then
        !	if (jstart.eq.-1) then
        	
        !	end if
!	        write(*,*) boundaryMeshName(k)
	        i=iconvname(boundaryMeshName(k))
!	        write(*,*) 'i',i
	        if (isdispls .ne. -1 ) then
	            sDispls(i+1)=sDispls(isdispls+1)+recvCounts(isdispls+1) 
!	            sDisplsp(i+1)=sDisplsp(isdispls+1)+recvCountsp(isdispls+1) 
	        end if
	        isdispls=i
	        recvCounts(i+1)=boundaryMeshSize(k)
!	        recvCountsp(i+1)=nodeboundaryMeshSize(k)*9
	    end if
	end do
	
!	do i = 1,nump
!	    write(*,*) '@@@@ rank=',rank, ' recvCounts=',recvCounts(i),' sDispls=',sDispls(i)
!	end do
	
		
!	m = recvCounts(1)
!	do i = 2,nump
!	    if(recvCounts(i).ne.0) then
!	        sDispls(i) = m
!	        m = m + recvCounts(i)
!	    end if
!	end do

!	write(*,*) 'rank,jpropeller=',rank,jpropeller
	do i=1,nump
        sendCounts(i)=recvCounts(i)
        rDispls(i)=sDispls(i)
!        sendCountsp(i)=recvCountsp(i)
!        rDisplsp(i)=sDisplsp(i)
            
!        sendCountsv(i)=15*sendCounts(i) ! ? 15 - leght of internal massive
!        recvCountsv(i)=15*recvCounts(i)
!        rDisplsv(i)=15*sDispls(i)
!        sDisplsv(i)=15*sDispls(i)

        sendCountsv(i)=18*sendCounts(i) ! ? 15 - leght of internal massive
        recvCountsv(i)=18*recvCounts(i)
        rDisplsv(i)=18*sDispls(i)
        sDisplsv(i)=18*sDispls(i)
	end do

!	write(*,*) '-------------------'
!	do i = 1,nump
!	    write(*,'(a5,i3,a12,i3,a2,i5,a8,i5)') 'proc=',rank, ' recvCounts[',i,']=',recvCounts(i), ' sDispls=',sDispls(i)
!	end do
!	write(*,*) '-------------------'


	netr0=0
!	netp0=0
	do i=1,nump
	    netr0=netr0+recvCounts(i)
!	    netp0=netp0+recvCountsp(i)
	end do

!	allocate (send(15*netr0),recv(15*netr0))
	allocate (send(18*netr0),recv(18*netr0))
	
!	allocate (sendp(20*netp0),recvp(20*netp0))
!	write(*,*) 'rank,20*netp0=',rank,20*netp0
	allocate (opposingOwnerFaceLabel(nFaces))
	allocate (meshCells(6,nCells))
	! end. MPI + rotating coor system (propeller)

	! set connectivities
	meshCells=0	
	do i=1,nFaces
	    do j=1,6
	        if (meshCells(j,meshFaceOwner(i)).eq.0) then
	            meshCells(j,meshFaceOwner(i))=i	
	            exit ! goto 1
	        end if
	    end do
! 1	    continue
	end do
	
	do i=1,nInternalFaces
	    do j=1,6
	        if (meshCells(j,meshFaceNeighbour(i)).eq.0) then
	            meshCells(j,meshFaceNeighbour(i))=i	
	            exit ! goto 2
	        end if
	    end do
! 2	    continue
	end do	
	
	do i=1,nFaces
	    do j=1,6
	        m=0
	        if (meshCells(j,meshFaceOwner(i)).ne.i) then
	            do k=1,4
	                do l=1,4
	                    if (meshFaces(k,i) &
                            .eq.meshFaces(l,meshCells(j,meshFaceOwner(i) ) ) ) m=1	
	                end do
	            end do	
	            if (m.eq.0) then
	                opposingOwnerFaceLabel(i)=meshCells(j,meshFaceOwner(i))
	                goto 3
	            end if
	        end if
	    end do
 3	    continue
        if (opposingOwnerFaceLabel(i).eq.0)then
            write(*,*) '@Fatal error: Face=',i
            stop 'Fatal error in opposingOwnerFaceLabel(i)'
        end if
	end do
	
	allocate (opposingNeighbourFaceLabel(nInternalFaces))
	do i=1,nInternalFaces
	    do j=1,6
	        m=0
	        if (meshCells(j,meshFaceNeighbour(i)).ne.i) then
	            do k=1,4
	                do l=1,4
	                    if (meshFaces(k,i) &
                            .eq.meshFaces(l,meshCells(j,meshFaceNeighbour(i) ) ) ) m=1	
	                end do
	            end do	
	            if (m.eq.0) then
	                opposingNeighbourFaceLabel(i)=meshCells(j,meshFaceNeighbour(i))
	                goto 4
	            end if
	        end if
	    end do
 4	    continue
	end do	
	! end of set connectivities

	! geometry setup
	allocate (meshCellPoints(8,nCells))
	meshCellPoints=0
	do i=1,nFaces
	    j=meshFaceOwner(i)
        call add_1(meshCellPoints(1:8,j),meshFaces(1:4,i))
    end do
	do i=1,nInternalFaces
	    j=meshFaceNeighbour(i)
        call add_1(meshCellPoints(1:8,j),meshFaces(1:4,i))
    end do
	open(32,file=processor(1:ie)//'outmesh',form='unformatted')
	write(32) nPoints,nCells,((meshPoints(jj,i),jj=1,3),i=1,nPoints) &
      ,((meshCellPoints(jj,i),jj=1,8),i=1,nCells)
    close(32)
    
	allocate (meshSf(3,nFaces))

!    open(2,file=processor(1:ie)//'meshSf')
!    write(2,*) ' i Fx Fy Fz x y z'
	do i=1,nFaces
	    call sxdef(meshSf(1:3,i),meshPoints(1:3,meshFaces(1,i)) & ! calculate normal arias
	                        ,meshPoints(1:3,meshFaces(2,i)) &
	                        ,meshPoints(1:3,meshFaces(3,i)) &
	                        ,meshPoints(1:3,meshFaces(4,i)) )

!	    write(2,'(i3,3e13.5,3f8.3)') i, meshSf(1:3,i), xg(1:3)
	    
	    
    !	write(*,*) 'Face,i',i,meshSf(1:3,i)
	    xg(1:3)=0.0d0
	    l=meshFaceOwner(i)
	    ! calc. xg - vector from center of cell to face
	    do j=1,8
	        xg(1:3)=xg(1:3)-meshPoints(1:3,meshCellPoints(j,l))/8.0d0
	    end do
	    do j=1,4
	        xg(1:3)=xg(1:3)+meshPoints(1:3,meshFaces(j,i))/4.0d0
	    end do
	    r=xg(1)*meshSf(1,i)+xg(2)*meshSf(2,i)+xg(3)*meshSf(3,i) ! projection of xg on normal of face
	    if (r.lt.0.0d0) meshSf(1:3,i)=-meshSf(1:3,i)

!	    xg(1:3) = 0.25 * (meshPoints(1:3,meshFaces(1,i)) + meshPoints(1:3,meshFaces(2,i)) &
!	                        + meshPoints(1:3,meshFaces(3,i)) &
!	                        + meshPoints(1:3,meshFaces(4,i)))

!	    write(2,'(i3,3e13.5,3f8.3)') i, meshSf(1:3,i), xg(1:3)
	end do
!	close(2)

! insert for check 20101203
!	meshSf=-meshSf
    ! set meshV, meshFaceCenters, meshCellCenters
	allocate (meshV(nCells),meshFaceCenters(3,nFaces),meshCellCenters(3,nCells))
	meshCellCenters=0.0d0
	meshV=0.0d0
	do i=1,nCells
	    xg(1:3)=0.0d0
	    do j=1,8
	        xg(1:3)=xg(1:3)-meshPoints(1:3,meshCellPoints(j,i))/8.0d0
	    end do
	    meshCellCenters(1:3,i)=xg(1:3)	
	end do
	do k=1,nFaces
	    i=meshFaceOwner(k)	
	    meshV(i)=meshV(i)+meshSf(1,k)                      &
	                           *(meshPoints(1,meshFaces(1,k)) &
	                            +meshPoints(1,meshFaces(2,k)) &
	                            +meshPoints(1,meshFaces(3,k)) &
	                            +meshPoints(1,meshFaces(4,k)) )
	    meshFaceCenters(1:3,k)= (meshPoints(1:3,meshFaces(1,k)) &
	                            +meshPoints(1:3,meshFaces(2,k)) &
	                            +meshPoints(1:3,meshFaces(3,k)) &
	                            +meshPoints(1:3,meshFaces(4,k)) )/4.0d0
	end do
	do k=1,nInternalFaces
	    i=meshFaceNeighbour(k)	
	    meshV(i)=meshV(i)-meshSf(1,k)                      &
	                       *(meshPoints(1,meshFaces(1,k)) &
	                        +meshPoints(1,meshFaces(2,k)) &
	                        +meshPoints(1,meshFaces(3,k)) &
	                        +meshPoints(1,meshFaces(4,k)) )
	end do
	vol=0.0d0
	Vmax=-1d10
	do i=1,nCells
    	meshV(i)=dabs(meshV(i)/4.0d0)
	    if (Vmax.lt.dabs(meshV(i))) Vmax=dabs(meshV(i))
	    vol=vol+meshV(i)
!	write(*,*) 'i=',i,meshV(i)
	end do
	
	Vmin = minVal(meshV(1:nCells))

    call OneValueMinMax(Vmin,Vmax)
    if(rank.eq.0)then
    	write(*,*) 'Vmax,Vmin=',Vmax,Vmin ! write(*,*) 'rank,vol,Vmax,Vmin=',rank,vol,Vmax,Vmin
    end if
	
	if(rank.eq.0)then
	    write(*,*) 'Volumes have been computed successfully!'
	end if

	allocate (Db(nFaces))
	allocate (Df(nFaces))
	Db=0.0d0
	Df=0.0d0
	Dbmin=1d10
	Dfmin=1d10
	do i=1,nFaces
	    j=opposingOwnerFaceLabel(i)
	    do k=1,3
	        Db(i)=Db(i)+(              &
                  (meshPoints(k,meshFaces(1,i)) &
                 +meshPoints(k,meshFaces(2,i)) &
                 +meshPoints(k,meshFaces(3,i)) &
                 +meshPoints(k,meshFaces(4,i)))/4.0d0 &
                 -(meshPoints(k,meshFaces(1,j)) &
                 +meshPoints(k,meshFaces(2,j)) &
                 +meshPoints(k,meshFaces(3,j)) &
                 +meshPoints(k,meshFaces(4,j)))/4.0d0 )**2
	    end do
	    Db(i)= dsqrt(Db(i))
	    if (Db(i).lt.Dbmin) Dbmin=Db(i)
	end do	
 !   write(*,*) 'rank,Dbmin=',rank,Dbmin

	if(rank.eq.0)then
	    write(*,*) 'Db-s have been computed successfully!'
	end if

	do i=1,nInternalFaces
	    j=opposingNeighbourFaceLabel(i)
	    do k=1,3
	        Df(i)=Df(i)+(              &
                  (meshPoints(k,meshFaces(1,i)) &
                 +meshPoints(k,meshFaces(2,i)) &
                 +meshPoints(k,meshFaces(3,i)) &
                 +meshPoints(k,meshFaces(4,i)))/4.0d0 &
                 -(meshPoints(k,meshFaces(1,j)) &
                 +meshPoints(k,meshFaces(2,j)) &
                 +meshPoints(k,meshFaces(3,j)) &
                 +meshPoints(k,meshFaces(4,j)))/4.0d0 )**2
	    end do
	    Df(i)= dsqrt(Df(i))
	    if (Df(i).lt.Dfmin) Dfmin=Df(i)	
	end do

	if(rank.eq.0)then
	    write(*,*) 'Df-s have been computed successfully!'
	end if

    ! end of geometry setup
!    write(*,*) 'rank,Dfmin=',rank,Dfmin

!	Df(nInternalFaces+1:nFaces)=Db(nInternalFaces+1:nFaces)
    ! MPI, geometru setup
!	m=0
!	do k=1,nBoundaryFields
!	    if (boundaryMeshName(k)(1:4).eq.'proc') then
!	        do i=1,boundaryMeshSize(k)
!	            m=m+1
!	            send(m)=Db(i+boundaryMeshStart(k))
!	        end do
!	    end if

!	end do
!    call MPI_Alltoallv(send,sendCounts,sDispls,MPI_DOUBLE_PRECISION,&
!                       recv,recvCounts,rDispls,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,ier)
!	m=0
!	do k=1,nBoundaryFields
!	    if (boundaryMeshName(k)(1:4).eq.'proc') then
!	        do i=1,boundaryMeshSize(k)
!	            m=m+1
!	            Df(i+boundaryMeshStart(k))=recv(m)
!	        end do
!	    end if
!	end do

    do k=1,nBoundaryFields
     if (boundaryMeshName(k)(1:4).eq.'proc') then
       iProc=iconvname(boundaryMeshName(k)) ! rank of neighbour proc
       call MPI_SENDRECV(Db(1+boundaryMeshStart(k)), boundaryMeshSize(k), &
                                        MPI_DOUBLE_PRECISION, iProc, 0, &
                         Df(1+boundaryMeshStart(k)), boundaryMeshSize(k), &
                                        MPI_DOUBLE_PRECISION, iProc, 0, &
                                         MPI_COMM_WORLD, statusp, ier)
      end if
    end do
	
	
    allocate (meshVPoints(nPoints))
    meshVPoints=0.0d0
    do i=1,nCells
        do j=1,8
            k=meshCellPoints(j,i)
            meshVPoints(k)=meshVPoints(k)+meshV(i)/2.0d0
        end do
    end do
    
    
! periodic
	do j=1,nBoundaryFields
        if (btype(j).eq.5) then
	        ndeltap=boundaryMeshSize(j)/2
 	        do k=1+boundaryMeshStart(j),boundaryMeshStart(j)+ndeltap
	            i=meshFaceOwner(k)
	            j1=meshFaceOwner(k+ndeltap)
	            do l=1,4
            !	write(*,'(i5,100f12.5)') l,(meshPoints(jj,meshFaces(l,k)),jj=1,3)
            !	write(*,'(i5,100f12.5)') l,(meshPoints(jj,meshFaces(l,k+ndeltap)),jj=1,3)
	                meshVPoints(meshFaces(l,k))=meshVPoints(meshFaces(l,k)) + meshV(j1)/2.0d0
	                meshVPoints(meshFaces(l,k+ndeltap))=meshVPoints(meshFaces(l,k+ndeltap)) + meshV(i)/2.0d0
	            end do
	        end do
	    end if
	end do
! processor
!	include 'netInterfaceInit.H'
    do i=1,nodeInterface
       sendVPoints(i)=meshVPoints(nodeBoundaryMeshlabel(i))
    end do

    recvVPoints = 0

    do iProc = 0,nump-1

    IF (sendCountsVPoints(iProc+1)/=0)THEN
!! recCountsVPoints=sendCountsVPoints, sDisplsVPoints=rDisplsVPoints
!! could replace sDisplsVPoints with cumulative total from sendCountsVPoints??? 
    call MPI_SENDRECV(sendVpoints(1+sDisplsVPoints(iProc+1)), sendCountsVPoints(iProc+1), &
                                                          MPI_DOUBLE_PRECISION, iProc, 0, &
                      recvVPoints(1+sDisplsVPoints(iProc+1)), sendCountsVPoints(iProc+1), &
                                                          MPI_DOUBLE_PRECISION, iProc, 0, &
                                                             MPI_COMM_WORLD,statusp, ier)
    ENDIF
    end do

!    call MPI_Alltoallv(sendVPoints,sendCountsVPoints,sDisplsVPoints,MPI_DOUBLE_PRECISION &
!	        ,recvVPoints,sendCountsVPoints,sDisplsVPoints,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,ier)


    do i=1,nodeInterface
       meshVPoints(nodeBoundaryMeshlabel(i))=meshVPoints(nodeBoundaryMeshlabel(i))+ recvVPoints(i)
    end do
    ! end MPI, geometru setup
    
  	! initialiozation DbCells - length of cells
	allocate (DbCells(nCells))
	do i=1,nCells
	    ds1=1d10
	    do j=1,6
    	    !if (Db(meshCells(j,i)).lt.ds1) ds1=Db(meshCells(j,i))
    	    iFace = meshCells(j,i)
    	    if(i.ne.meshFaceOwner(iFace)) cycle
    	    ds1 = min(ds1,Db(iFace))
!    	    if(iFace.le.nInternalFaces) then
!    	        ds1 = min(ds1,Df(iFace))
!    	    end if
	    end do
	    DbCells(i)=ds1
	end do
	
	allocate(Points2Faces(FacesInPointsMax,nPoints))
	Points2Faces = 0
	
	do iFace = 1,nFaces
	    do PointNum = 1,4
	        iPoint = meshFaces(PointNum,iFace)
	        do FaceNum = 1,FacesInPointsMax
	            if(Points2Faces(FaceNum,iPoint).eq.0)then
	                Points2Faces(FaceNum,iPoint) = iFace
	                exit
	            end if
	        end do
	    end do
	end do
	
	! checking cyclic boundaries
	call CheckCyclicBC()

	if(rank.eq.0)then
	    write(*,*) 'AfterReadMeshOpenFOAM has been started ...'
	end if
	
	call AfterReadMeshOpenFOAM()

	if(rank.eq.0)then
	    write(*,*) 'AfterReadMeshOpenFOAM has been finished successfully!'
	end if
	
end subroutine ReadMeshOpenFOAM

subroutine CheckCyclicBC()
use Data_Cfoam, only: btype
use Data_Mesh
implicit none

    integer :: j, iFace,jFace, bmeshs, ndeltap,bmeshf
    real(8) :: iR(3),jR(3), rr

    do j=1,nBoundaryFields
        if (btype(j).ne.5) cycle ! goto 10
        bmeshs = 1 + boundaryMeshStart(j) ! Phil, 2012-02-28
        ndeltap = boundaryMeshSize(j)/2
        bmeshf = boundaryMeshStart(j) + ndeltap
        
        do iFace = bmeshs,bmeshf
            jFace = iFace+ndeltap ! opposingOwnerFaceLabel(iFace+ndeltap)
            iR = 0d0
            iR(1:2) = meshFaceCenters(1:2,iFace)
            jR = 0d0
            jR(1:2) = meshFaceCenters(1:2,jFace)
            iR = iR - jR
            
            rr = sqrt(sum(iR*iR))

            if(rr.gt.1d-5)then
                write(*,*) 'CyclicBC is worng:'
                write(*,*) 'iFace =', iFace
                write(*,'(a4,3f12.5)') 'iR =', iR
                write(*,*) 'jFace =', jFace
                write(*,'(a4,3f12.5)') 'jR =', jR
            end if
            
            iR(1:3) = meshsf(1:3,iFace)
            jR(1:3) = meshsf(1:3,jFace)
            
            rr = abs(dsqrt(sum(iR*iR)) - dsqrt(sum(jR*jR)))

            if(rr.gt.1d-10)then
                write(*,*) 'CyclicBC areas are worng:'
                write(*,*) 'iFace =', iFace
                write(*,'(a4,3f12.5)') 'iS =', dsqrt(sum(iR*iR))
                write(*,*) 'jFace =', jFace
                write(*,'(a4,3f12.5)') 'jS =', dsqrt(sum(jR*jR))
            end if
        end do

    end do
    
!    stop 'CheckCyclicBC has done!'

end subroutine CheckCyclicBC


subroutine CalcNewPeriodicNumeration()
use Data_cfoam, only: PeriodicBound1, PeriodicBound2
use Data_Mesh
use Data_cfoam, only: ZeroEpsilon
implicit none
! assumed that PeriodicBound2 = PeriodicBound1 + 1, and these last boundaries in boundary
! moreover both boundaries are oppostite sides of 1 layer cells

    integer :: iCell,iFace, jFace, nPeriodic, i,j, iFaceStart,iFaceEnd, iFace1, i2j(4), iPoint,jPoint
    integer,allocatable :: iFaceNew2iFace(:), iFace2iFaceNew(:)
    integer,allocatable :: AdjacentFace(:) ! for periodic boundary
	integer,allocatable :: boundaryMeshStartNew(:)

    integer :: nPeriodicFaces,nPeriodicFaces2, bmeshs,bmeshf
    integer :: iFaceGlobal, iAdjacentFace, iAdjacentFaceGlobal
    real(8) :: FaceCenter(3), FaceDistance(3), Pc(3)

    nPeriodicFaces = boundaryMeshSize(PeriodicBound1)
    nPeriodicFaces2 = boundaryMeshSize(PeriodicBound2)

    if(nPeriodicFaces.ne.nPeriodicFaces2.or.PeriodicBound1.gt.PeriodicBound2)then
        Stop '@@ Fatal Error: if(nFacesPeriodic1.ne.nFacesPeriodic2.or.PeriodicBound.gt.PeriodicBound2) !'
    end if

!    allocate(iFaceNew2iFace(nFaces),iFace2iFaceNew(nFaces))
    allocate(iFace2iFaceNew(nPeriodicFaces))

    bmeshs = 1 + boundaryMeshStart(PeriodicBound1) ! Phil, 2012-02-28
    bmeshf = boundaryMeshStart(PeriodicBound1) + nPeriodicFaces

    i2j(1) = 1
    i2j(2) = 4
    i2j(3) = 3
    i2j(4) = 2

!    do iFace = bmeshs,bmeshf
    do iFace = 1,nPeriodicFaces
        iFace1 = iFace+boundaryMeshStart(PeriodicBound1)
        jFace = opposingOwnerFaceLabel(iFace+boundaryMeshStart(PeriodicBound1))
        
        do j = 1,4
            iPoint = meshFaces(j,iFace1)
            jPoint = meshFaces(i2j(j),jFace)
            Pc(1:2) = 0.5*(meshPoints(1:2,iPoint) + meshPoints(1:2,jPoint))
            
            meshPoints(1:2,iPoint) = Pc(1:2)
            meshPoints(1:2,jPoint) = Pc(1:2)
        end do

        iFace2iFaceNew(iFace) = jFace
    end do

!!!    allocate(AdjacentFace(nPeriodicFaces))
!!!    allocate(boundaryMeshStartNew(nBoundaryFields))
!!    
!!    do iFace = 1,nFaces
!!        iFaceNew2iFace(iFace) = iFace
!!    end do
!!    
!!    ! Global shift
!!    boundaryMeshStartNew = boundaryMeshStart
!!    do j = PeriodicBound1+1, PeriodicBound2-1
!!        iFaceStart = boundaryMeshStartNew(j) + 1
!!        iFaceEnd = boundaryMeshStartNew(j) + boundaryMeshSize(j)
!!        iFaceNew2iFace(iFaceStart:iFaceEnd) = iFaceNew2iFace(iFaceStart:iFaceEnd) - nPeriodicFaces
!!
!!        boundaryMeshStartNew(j) = boundaryMeshStart(j) - nPeriodicFaces
!!    end do
!!
!!    iFaceStart = boundaryMeshStart(PeriodicBound1) + 1
!!    iFaceEnd = boundaryMeshStart(PeriodicBound1) + boundaryMeshSize(PeriodicBound1)
!!    iFaceNew2iFace(iFaceStart:iFaceEnd) = iFaceNew2iFace(iFaceStart:iFaceEnd) + &
!!                                          nFaces - 2*nPeriodicFaces - boundaryMeshStartNew(PeriodicBound1)
!!
!!    boundaryMeshStartNew(PeriodicBound1) = nFaces - 2*nPeriodicFaces
!!
!!    do j = PeriodicBound2+1, nBoundaryFields
!!        iFaceStart = boundaryMeshStartNew(j) + 1
!!        iFaceEnd = boundaryMeshStartNew(j) + boundaryMeshSize(j)
!!        iFaceNew2iFace(iFaceStart:iFaceEnd) = iFaceNew2iFace(iFaceStart:iFaceEnd) - 2*nPeriodicFaces
!!
!!        boundaryMeshStartNew(j) = boundaryMeshStart(j) - 2*nPeriodicFaces
!!    end do
!!
!!    iFaceStart = boundaryMeshStart(PeriodicBound2) + 1
!!    iFaceEnd = boundaryMeshStart(PeriodicBound2) + boundaryMeshSize(PeriodicBound2)
!!    iFaceNew2iFace(iFaceStart:iFaceEnd) = iFaceNew2iFace(iFaceStart:iFaceEnd) + &
!!                                          nFaces - nPeriodicFaces - boundaryMeshStartNew(PeriodicBound2)
!!
!!    boundaryMeshStartNew(PeriodicBound2) = nFaces - nPeriodicFaces
!!    ! End of Global shift, now Periodic boundaries are lie in the end of array.
!!
!!    ! Definition adjacent faces for periodic boundaries
!!    AdjacentFace = 0
!!    do iFace = 1,nPeriodicFaces
!!        iFaceGlobal = boundaryMeshStart(PeriodicBound1) + iFace
!!        FaceCenter(1:3) = meshFaceCenters(1:3,iFaceGlobal)
!!        do iAdjacentFace = 1,nPeriodicFaces
!!            iAdjacentFaceGlobal = boundaryMeshStart(PeriodicBound2) + iAdjacentFace
!!            FaceDistance(1:3) = meshFaceCenters(1:3,iAdjacentFaceGlobal) - FaceCenter
!!            
!!            if(sum(FaceDistance*FaceDistance).lt.ZeroEpsilon)then
!!                AdjacentFace(iFace) = iAdjacentFace
!!            end if
!!        end do
!!    end do
!!    
!!    do iFace = 1,nPeriodicFaces
!!        write(1,*) iFace, AdjacentFace(iFace)
!!    end do
!!
!!    iFaceStart = boundaryMeshStart(PeriodicBound2) + 1
!!    iFaceEnd = boundaryMeshStart(PeriodicBound2) + boundaryMeshSize(PeriodicBound2)
!!    iFaceNew2iFace(iFaceStart:iFaceEnd) = nFaces - nPeriodicFaces + AdjacentFace(1:nPeriodicFaces)
!!    
!!    do iFace = 1,nFaces
!!        iFace2iFaceNew(iFaceNew2iFace(iFace)) = iFace
!!    end do
!!    ! End of Definition adjacent faces for periodic boundaries
    
    ! output new numeration
    
    boundaryMeshSize(PeriodicBound1) = 2 * boundaryMeshSize(PeriodicBound1)
    
!    ! boundary_new
!        open(2,file='boundary')
!        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
!        write(2,*) '| =========                 |                                                 |'
!        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
!        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
!        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
!        write(2,*) '|    \\/     M anipulation  |                                                 |'
!        write(2,*) ' ' ! \*---------------------------------------------------------------------------*/'
!        write(2,*) 'FoamFile'
!        write(2,*) '{'
!        write(2,*) '    version     2.0;'
!        write(2,*) '    format      ascii;'
!        write(2,*) '    class       polyBoundaryMesh;'
!        write(2,*) '    location    "constant/polyMesh";'
!        write(2,*) '    object      boundary;'
!        write(2,*) '}'
!        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
!        write(2,*) ' '
!        
!        write(2,'(i3)') nBoundaryFields-1
!        write(2,*) '('
!        
!        do i = 1,nBoundaryFields
!            if(i.eq.PeriodicBound1.or.i.eq.PeriodicBound2) cycle
!            call writeBoundaryProps(2,i)
!        end do
!        call writeBoundaryProps(2,PeriodicBound1)
!        
!        write(2,*) ')'
!        write(2,*) ' '
!        write(2,*) '// ************************************************************************* //'
!
!        close(2)

    ! faces_new
        open(2,file='faces')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) ' ' ! \*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       faceList;'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      faces;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) nFaces
        write(2,*) '('
        
        do iFace = 1,bmeshf ! nFaces
            write(2,*) '4(',meshFaces(1:4,iFace)-1,')'
        end do
        
        do iFace = 1,nPeriodicFaces
            write(2,*) '4(',meshFaces(1:4,iFace2iFaceNew(iFace))-1,')'
        end do
        
        write(2,*) ')'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'

        close(2)

    ! owner_new
        open(2,file='owner')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) ' ' ! \*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       labelList;'
        write(2,*) '    note        "nPoints:',nPoints, ' nCells:', nCells, ' nFaces:', &
                    nFaces, ' nInternalFaces:',nInternalFaces, '";'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      owner;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) nFaces
        write(2,*) '('

        do iFace = 1,bmeshf ! nFaces
            write(2,*) meshFaceOwner(iFace)-1
        end do
        
        do iFace = 1,nPeriodicFaces
            write(2,*) meshFaceOwner(iFace2iFaceNew(iFace))-1
        end do

        write(2,*) ')'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'

        close(2)

    ! points_new
        open(2,file='points')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) ' ' ! \*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       vectorField;'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      points;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        
        write(2,*) nPoints
        write(2,*) '('
        
        do iPoint = 1,nPoints
            write(2,*) '(',meshPoints(1:3,iPoint),')'
        end do
        
        write(2,*) ')'
        write(2,*) ' '
        write(2,*) '// ************************************************************************* //'

        close(2)
    ! end of output new numeration

    Stop 'NewPeriodicNumeration done!'

contains

    subroutine writeBoundaryProps(FileNum,i)
    implicit none
        integer,intent(in) :: FileNum,i

        write(2,*) boundaryMeshName(i)
        write(2,*) '    {'
        write(2,*) '        type ',boundaryMeshType(i),';'
        write(2,*) '        nFaces ', boundaryMeshSize(i), ';'
        write(2,*) '        startFace ', boundaryMeshStartNew(i), ';'
        write(2,*) '    }'

    end subroutine writeBoundaryProps

end subroutine CalcNewPeriodicNumeration

! assume that periodic boundary has number nBoundField
subroutine PR_SplitPeriodicBoundaries()
use Data_Mesh
use Data_cfoam, only: ZeroEpsilon
implicit none

    integer :: iCell,iFace, nPeriodic, i,j, iFaceStart,iFaceEnd
    integer :: nPeriodicFaces, nFacesNew, nInternalFacesNew, nPointsNew
    integer :: iFacePoint,jFacePoint, iPoint,jPoint
    logical :: PointFind
    integer :: bmeshs, bmeshf, iBound
	integer,allocatable :: boundaryMeshStartNew(:)
    integer,allocatable :: iFaceNew2iFace(:), iFace2iFaceNew(:)
    integer,allocatable :: PointIndicator(:),PointOrderNew(:),PointIndexNew(:)
    real(8) :: r(3), rMod

    allocate(PointIndicator(nPoints),PointOrderNew(nPoints),PointIndexNew(nPoints))

    nPeriodicFaces = boundaryMeshSize(nBoundaryFields)/2
    nFacesNew = nFaces-nPeriodicFaces
    nInternalFacesNew = nInternalFaces+nPeriodicFaces

    ! define new points
    do i = 1,nPoints
        PointIndicator(i) = i
        PointOrderNew(i) = i
    end do

    bmeshs = boundaryMeshStart(nBoundaryFields) + 1
    bmeshf = boundaryMeshStart(nBoundaryFields) + nPeriodicFaces
    do iFace = bmeshs,bmeshf
        do iFacePoint = 1,4
            iPoint = meshFaces(iFacePoint,iFace)
            PointFind = .false.
            do jFacePoint = 1,4
                jPoint = meshFaces(jFacePoint,iFace+nPeriodicFaces)
                
                r(1:3) = meshPoints(1:3,iFacePoint) - meshPoints(1:3,jFacePoint)
                rMod = sum(r(1:3)*r(1:3))
                if(rMod.lt.ZeroEpsilon)then
                    PointFind = .true.
                    PointIndicator(jPoint) = 0
                    PointOrderNew(jPoint) = iPoint
                    exit
                end if
            end do
            if(PointFind.eqv..false.)then
                write(*,*) '@@Fatala error adjacent points have not been found! iPoint=', &
                        iPoint, 'jPoint=',jPoint
                Stop '@@Fatal Error!'
            end if
        end do
    end do
    
    nPointsNew = 0
    do iPoint = 1,nPoints
        if(PointIndicator(iPoint).eq.0) cycle
        if(PointIndicator(iPoint).eq.-1)then
            write(*,*) '@@Fatala error PointIndicator(iPoint).eq.-1! iPoint=', iPoint
            Stop '@@Fatal Error!'
        end if
        nPointsNew = nPointsNew + 1
        PointIndexNew(iPoint) = nPointsNew
    end do

    ! end of definition new points

    allocate(boundaryMeshStartNew(nBoundaryFields-1))
    do iBound = 1,nBoundaryFields-1
        boundaryMeshStartNew(iBound) = boundaryMeshStart(iBound) + nPeriodicFaces
    end do
    
    allocate(iFaceNew2iFace(nFacesNew))
    
    do iFace = 1,nInternalFaces
        iFaceNew2iFace(iFace) = iFace
    end do
    do iFace = 1,nPeriodicFaces
        i = iFace + nInternalFaces
        iFaceNew2iFace(i) = boundaryMeshStart(nBoundaryFields) + iFace
    end do
    do iFace = nInternalFacesNew+1,nFacesNew
        iFaceNew2iFace(iFace) = iFace - nPeriodicFaces
    end do
    
    allocate(iFace2iFaceNew(nFacesNew))
    do iFace = 1,nFacesNew
        iFace2iFaceNew(iFace) = iFaceNew2iFace(iFace)
    end do

    ! boundary_new
        open(2,file='boundary')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) ' ' ! \*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       polyBoundaryMesh;'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      boundary;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        
        write(2,'(i3)') nBoundaryFields-1
        write(2,*) '('
        
        do i = 1,nBoundaryFields-1
            call writeBoundaryProps(2,i)
        end do
        
        write(2,*) ')'
        write(2,*) ' '
        write(2,*) '// ************************************************************************* //'

        close(2)

    ! faces_new
        open(2,file='faces')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) '\*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       faceList;'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      faces;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) nFacesNew
        write(2,*) '('
        
        do iFace = 1,nFacesNew
!            write(2,*) '4(',meshFaces(1:4,iFace2iFaceNew(iFace))-1,')'
            write(2,*) '4(',PointIndexNew(PointOrderNew(meshFaces(1:4,iFaceNew2iFace(iFace))))-1,')'
        end do
        
        write(2,*) ')'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'

        close(2)

    ! owner_new
        open(2,file='owner')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) ' ' ! \*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       labelList;'
        write(2,*) '    note        "nPoints:',nPointsNew, ' nCells:', nCells, &
                    ' nFaces:',nFacesNew, ' nInternalFaces:',nInternalFacesNew, '";'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      owner;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) nFacesNew
        write(2,*) '('

        do iFace = 1,nFacesNew
!            write(2,*) meshFaceOwner(iFace2iFaceNew(iFace))-1
            write(2,*) meshFaceOwner(iFaceNew2iFace(iFace))-1
        end do

        write(2,*) ')'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'

        close(2)


    ! neighbour_new
        open(2,file='neighbour')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) '\*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       labelList;'
        write(2,*) '    note        "nPoints:',nPointsNew, ' nCells:', nCells, &
                    ' nFaces:',nFacesNew, ' nInternalFaces:',nInternalFacesNew, '";'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      neighbour;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) nInternalFacesNew
        write(2,*) '('

!        do iFace = 1,nInternalFacesNew
!!            write(2,*) meshFaceOwner(iFace2iFaceNew(iFace))-1
!            write(2,*) meshFaceOwner(iFaceNew2iFace(iFace))-1
!        end do

        do iFace = 1,nInternalFaces
!            write(2,*) meshFaceNeighbour(iFace2iFaceNew(iFace))-1
!            write(2,*) meshFaceNeighbour(iFaceNew2iFace(iFace))-1
            write(2,*) meshFaceNeighbour(iFace)-1
        end do
        bmeshs = boundaryMeshStart(nBoundaryFields) + nPeriodicFaces + 1
        bmeshf = boundaryMeshStart(nBoundaryFields) + boundaryMeshSize(nBoundaryFields)
        do iFace = bmeshs,bmeshf
            write(2,*) meshFaceOwner(iFace)-1
        end do
        write(2,*) ')'
        write(2,*) ' '
        write(2,*) ' '
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'

        close(2)
    ! end of output new numeration

    ! points_new
        open(2,file='points')
        write(2,*) ' ' ! /*--------------------------------*- C++ -*----------------------------------*\'
        write(2,*) '| =========                 |                                                 |'
        write(2,*) '| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'
        write(2,*) '|  \\    /   O peration     | Version:  1.7.1                                 |'
        write(2,*) '|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |'
        write(2,*) '|    \\/     M anipulation  |                                                 |'
        write(2,*) ' ' ! \*---------------------------------------------------------------------------*/'
        write(2,*) 'FoamFile'
        write(2,*) '{'
        write(2,*) '    version     2.0;'
        write(2,*) '    format      ascii;'
        write(2,*) '    class       polyBoundaryMesh;'
        write(2,*) '    location    "constant/polyMesh";'
        write(2,*) '    object      points;'
        write(2,*) '}'
        write(2,*) '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
        write(2,*) ' '
        write(2,*) ' '
        
        write(2,*) nPointsNew
        write(2,*) '('
        
        do iPoint = 1,nPoints
            if(PointIndicator(iPoint).eq.0) cycle
            write(2,*) '(',meshPoints(1:3,PointOrderNew(iPoint)),')'
        end do
        
        write(2,*) ')'
        write(2,*) ' '
        write(2,*) '// ************************************************************************* //'

        close(2)

    Stop 'SplitPeriodicBoundaries done!'

contains

    subroutine writeBoundaryProps(FileNum,i)
    implicit none
        integer,intent(in) :: FileNum,i

        write(2,*) boundaryMeshName(i)
        write(2,*) '    {'
        write(2,*) '        type ',boundaryMeshType(i),';'
        write(2,*) '        nFaces ', boundaryMeshSize(i), ';'
        write(2,*) '        startFace ', boundaryMeshStartNew(i), ';'
        write(2,*) '    }'

    end subroutine writeBoundaryProps

end subroutine PR_SplitPeriodicBoundaries

subroutine AfterReadMeshOpenFOAM()
use Data_Mesh !, only: 
implicit none

    integer :: iCell, jCell, iFace, iCellFace, iCellNeighbours, nCellNeighboursMax

    allocate(nCellNeighbours(nCells))
    
    nCellNeighbours(1:nCells) = 0

    do iCell = 1,nCells
        do iCellFace = 1,6
            iFace= meshCells(iCellFace,iCell)
            if(iFace.gt.nInternalFaces) cycle
            nCellNeighbours(iCell) = nCellNeighbours(iCell) + 1
        end do
    end do

    nCellNeighboursMax = maxVal(nCellNeighbours(1:nCells))
    allocate(CellNeighbours(nCellNeighboursMax,nCells))
    
    CellNeighbours(1:nCellNeighboursMax,1:nCells) = 0
    
    do iCell = 1,nCells
        iCellNeighbours = 0
        do iCellFace = 1,6
            iFace= meshCells(iCellFace,iCell)
            if(iFace.gt.nInternalFaces) cycle
            jCell = meshFaceOwner(iFace)
            if(jCell.eq.iCell)then
                jCell = meshFaceNeighbour(iFace)
            end if
            iCellNeighbours = iCellNeighbours + 1

            CellNeighbours(iCellNeighbours,iCell) = jCell
        end do
    end do
    

end subroutine AfterReadMeshOpenFOAM


	subroutine inputMain(name,nPoints,nCells,nFaces,nInternalFaces)
    implicit none
        integer i,j
        character (*) name
	    character *280 a
        integer nPoints,nCells,nFaces,nInternalFaces
        
	    open(11,file=name)
	    do i=1,12
	        read(11,*)
	    end do
	    read(11,'(a280)') a
	    do j=1,250
	        if (a(j:j).eq.'"') a(j:j)=' '
	    end do
	    do j=1,250
	        if (a(j:j+7).eq.'nPoints:') read(a(j+8:280),*) nPoints
	        if (a(j:j+6).eq.'nCells:') read(a(j+7:280),*) nCells
	        if (a(j:j+6).eq.'nFaces:') read(a(j+7:280),*) nFaces
	        if (a(j:j+14).eq.'nInternalFaces:') read(a(j+15:280),*) nInternalFaces
	    end do
	    
!	    write(*,*) 'nPoints:',nPoints
!	    write(*,*) 'nCells:',nCells
!	    write(*,*) 'nFaces:',nFaces
!	    write(*,*) 'nInternalFaces:',nInternalFaces
	    
	    close(11)
	return
	end subroutine inputMain
	
	subroutine inputBoundary(name,nBoundaryFields)
!	use Data_cfoam, only: rank
    implicit none
        character (*) name
        integer nBoundaryFields
        integer i,j
	
	    open(11,file=name)
	    do i=1,17
	        read(11,*)
	    end do
	    
	    read(11,*) nBoundaryFields
	    
	    ! write(*,*) 'nBoundaryFields:',nBoundaryFields !, ' rank=', rank
	    
	    close(11)
	    return
	end subroutine inputBoundary
	    
	subroutine inputBoundary1(name,nBoundaryFields,boundaryMeshSize,boundaryMeshStart,boundaryMeshName &
	,boundaryMeshType)
        implicit none
	
	    character (*),intent(in) :: name
	    integer,intent(in) :: nBoundaryFields
	    integer,intent(out) :: boundaryMeshSize(nBoundaryFields),boundaryMeshStart(nBoundaryFields)
	    character(LEN=50),intent(out) :: boundaryMeshName(nBoundaryFields)
	    character(LEN=20),intent(out) :: boundaryMeshType(nBoundaryFields)	    
	    character (LEN=280) a
        integer i,j
        
	    open(11,file=name)
	    do i=1,19
	        read(11,*)
	    end do
	    do i=1,nBoundaryFields
	        read(11,*) boundaryMeshName(i)
!	        write(*,*) 'boundaryMeshName(i)',i,boundaryMeshName(i)
	        read(11,*)
        !	read(11,*)
	        read(11,'(a280)') a
	        do j=1,250
	            if (a(j:j).eq.';') a(j:j)=' '
	            if (a(j:j+3).eq.'type') a(j:j+3)='    '
	        end do
	        read (a,*) boundaryMeshType(i)
    	
	        read(11,'(a280)') a
	        do j=1,250
    	        if (a(j:j).eq.';') a(j:j)=' '
	        end do
	        read (a(21:280),*) boundaryMeshSize(i)
!	        write(*,*) 'boundaryMeshSize(i)',i,boundaryMeshSize(i)
	        read(11,'(a280)') a
	        do j=1,250
	            if (a(j:j).eq.';') a(j:j)=' '
	        end do
	        read (a(21:280),*) boundaryMeshStart(i)
!	        write(*,*) 'boundaryMeshStart(i)',i,boundaryMeshStart(i)
	        !boundaryMeshStart(i)=boundaryMeshStart(i)+1	
!	        write(*,*) boundaryMeshName(i)(1:4)
	        if (boundaryMeshType(i)(1:4).eq.'cycl') read(11,*)

	        if (boundaryMeshName(i)(1:4).eq.'proc') then
	            read(11,*)
	            read(11,*)
	        end if
	        read(11,*)
	    end do
	    
	    close(11)
	    return
	end subroutine inputBoundary1
	
	subroutine inputMain1(name,nInternalFaces,meshFaceOwner)
    implicit none
        character (*) name
	    character *280 a
        integer nInternalFaces
	    integer meshFaceOwner(1:nInternalFaces)
        integer i,j
        
	    open(11,file=name)
	    do i=1,21
	        read(11,*)
	    end do
	    do i=1,nInternalFaces
	        read(11,*) meshFaceOwner(i)
	        meshFaceOwner(i)=meshFaceOwner(i)+1
	    end do
	    
	    close(11)
!	    write(*,*) 'read '//name//' was done',nInternalFaces
	    return
	end subroutine inputMain1
	
	subroutine inputMain2(name,nInternalFaces,meshFaceOwner)
    implicit none
        character (*) name
	    character *280 a
        integer nInternalFaces
	    integer meshFaceOwner(4,nInternalFaces)
        integer i,j,ii
        
	    open(11,file=name)
	    do i=1,20
	        read(11,*)
	    end do
	    do i=1,nInternalFaces
	        read(11,'(a280)') a
	        do j=1,250
	            if (a(j:j).eq.')') a(j:j)=' '
	            if (a(j:j).eq.'(') a(j:j)=' '
	        end do
	        read(a,*) ii,(meshFaceOwner(j,i),j=1,4)
	        do j=1,4
	            meshFaceOwner(j,i)=meshFaceOwner(j,i)+1
	        end do
	    end do
	    close(11)
!	    write(*,*) 'read faces '//name//' was done',nInternalFaces
	    return
	end subroutine inputMain2
	
	subroutine inputMain3(name,nInternalFaces,meshFaceOwner)
    implicit none
        character (*) name
    	character *280 a
        integer nInternalFaces
	    real *8 meshFaceOwner(3,nInternalFaces)
        integer i,j
	    open(11,file=name)

	    do i=1,20
    	    read(11,*)
	    end do
	    do i=1,nInternalFaces
	        read(11,'(a280)') a
	        do j=1,250
	            if (a(j:j).eq.')') a(j:j)=' '
	            if (a(j:j).eq.'(') a(j:j)=' '
	        end do
	        read(a,*) (meshFaceOwner(j,i),j=1,3)
	    end do
	    
	    close(11)
!	    write(*,*) 'read faces '//name//' was done',nInternalFaces
	    return
	end subroutine inputMain3

    integer function iconvname(strx)
    use Data_cfoam, only: rank
    implicit none
!	    integer iconvname
        integer :: i,j
	    character(*), intent(in) :: strx

!	    do i=1,20
        i = 0
        do
            i = i + 1
            if (strx(i:i+1).eq.'to') then
	            j=i+2
	            exit
	        else
	            if(strx(i:i).eq.' ') then
                    write(*,*) '@@Fatal Error: wrong boundaryName, proc ', rank
                    Stop '@@Fatal Error: wrong boundaryName'
                end if
	        end if
	    end do
!	write(*,*) 'j,strx=',j,strx
        i = LEN_TRIM(strx)
	    read(strx(j:i),*) iconvname
	    return
	end function iconvname

	subroutine sxdef(sx,x1,x2,x3,x4)
	implicit none
	! in: x1,x2,x3,x4 - verticies of face
	! out: sx - face normal, lenght = face area
    !	IMPLICIT REAL *8(a-h,o-z)
	    real(8),intent(out) :: sx(3)
	    real(8),intent(in) :: x1(3)
	    real(8),intent(in) :: x2(3)
	    real(8),intent(in) :: x3(3)
	    real(8),intent(in) :: x4(3)
	    integer rank
    !	SX=0.5*((Y3-Y1)*(Z4-Z2)-(Y4-Y2)*(Z3-Z1))
    !	SY=0.5*((Z3-Z1)*(X4-X2)-(Z4-Z2)*(X3-X1))
    !	SZ=0.5*((X3-X1)*(Y4-Y2)-(X4-X2)*(Y3-Y1))
	    sx(1)=0.5d0*((x3(2)-x1(2))*(x4(3)-x2(3))         &
         -(x4(2)-x2(2))*(x3(3)-x1(3)))
	    sx(2)=0.5d0*((x3(3)-x1(3))*(x4(1)-x2(1))-(x4(3) &
         -x2(3))*(x3(1)-x1(1)))
	    sx(3)=0.5d0*((x3(1)-x1(1))*(x4(2)-x2(2)) &
         -(x4(1)-x2(1))*(x3(2)-x1(2)))
	    return
	end subroutine sxdef
